load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            User Config                   ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


if (.not. isvar("areaMaskFile")) then
    areaMaskFile = "/public/users/wrf/Data/WuJiang_Map/wujiang_domaian_mask_d02.csv"
    print("areaMaskFile is undefined, using default: " + areaMaskFile)
else
    print("areaMaskFile is: " + areaMaskFile)
end if

if (.not. isvar("input_file")) then
    input_file = "/public/users/wrf/wrfout/2012111408/wrfout_d02_2012-11-14_00:00:00"
    print("input_file is undefined, using default: " + input_file)
end if

if (.not. isvar("output_prefix")) then
    print("output_prefix is undefined")
    output_prefix = "~/wrfpost/test/areal_precip"
    print("output_prefix is undefined, using default: " + output_prefix)
    system("mkdir " + output_prefix)
end if

; default precipitation delta-t is 1 hour
if (.not. isvar("precip_dt")) then
    print("precip_dt is undefined")
    precip_dt = 1 ; output 1h precip
end if


; Read areaMaskFile
region_name = readAsciiHead(areaMaskFile, 1)
wrf_mask = readAsciiTable(areaMaskFile,4,"float",1)
ij_masks = tointeger(wrf_mask(:,2:3))
dim_mask = dimsizes(ij_masks)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;           Open NetCDF File               ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
a = addfile(input_file + ".nc", "r")
filename = systemfunc("basename " + input_file)

domain = substring(filename, 7, 9)
start_date = substring(filename, 11, 23)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;              Process Time                ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
times = a->Times
hours = WRF_Times2Udunits_c(times, 0)
ntimes = dimsizes(hours)
wrf_dt = (hours(ntimes-1) - hours(0))/(ntimes-1)

; UTC -> CST
hours = hours + 8 ; 
timeStr_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")

dt = round(precip_dt / wrf_dt, 3)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            Read Precipitation            ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
rain_exp = wrf_user_getvar(a,"RAINNC",-1)
rain_con = wrf_user_getvar(a,"RAINC",-1)
rain_tot = rain_exp + rain_con
dim_rain = dimsizes(rain_tot)

rain_crop = new((/dim_rain(0),dim_mask(0)/), "float")
do i = 0,dim_mask(0)-1,1
    rain_crop(:,i) = rain_tot(:,ij_masks(i,1), ij_masks(i,0))
end do
rain_crop_avg = dim_sum(rain_crop) / tofloat(dim_mask(0))
rain_tend = array_append_record(0.0, (rain_crop_avg(dt:ntimes-1:dt) - rain_crop_avg(0:ntimes-1-dt:dt)), 0)
timeStr_tend = timeStr_cst(0:ntimes-1:dt)
printArray = region_name + ", " + timeStr_tend + ", " + sprintf("%6.3f", rain_tend) + ", " + sprintf("%7.3f", rain_crop_avg(0:ntimes-1:dt))
header = "Region Name, Time, Precip Tendency, Total Precip"
printArray_all = array_append_record(header, printArray, 0)

output_filename = output_prefix + "/areal_precip_" + region_name + "_" + domain + "_" + start_date + ".csv"

asciiwrite(output_filename, printArray_all)
